function [hs, lats, lons] = CoordPath(s, perigee)
%
% Numerical Integration with Gauss Algorithm
%
%DESCRIPTION:
%This function implements the numerical integration of the TEC with Slat
%path using Gauss Algorithm. [2.5.8.2.8]
%
%PROTOTYPE:
% [hs, lats, lons] = CoordPath(s, perigee)
%
%--------------------------------------------------------------------------
% INPUTS:
%   s         [1x1]       Distance                  [km]
%   perigee   [---]       Perigee Results           [strc] (see NOTES)
%--------------------------------------------------------------------------
% OUTPUTS:
%   hs        [1x1]       Point P Height            [km]
%   lats      [1x1]       Point P Latitude          [deg]
%   lons      [1x1]       Point P Longitude         [deg]
%--------------------------------------------------------------------------
%
%NOTES:
% - The input "perigee" has been chosen to be a structure (for compactness of
%   the code) defined as:
%       perigee.rp      = Radius           [km]
%       perigee.latp    = Latitude         [deg]
%       perigee.lonp    = Longitude        [deg]
%       perigee.sinlatp = Sine of Latitude [-]
%       perigee.coslatp = Cosine of Latit. [-]
%       perigee.sinsigp = Sine of Zenith   [-]
%       perigee.cossigp = Cosine of Zenith [-]
%
%CALLED FUNCTIONS:
% (none)
%
%UPDATES:
% (none)
%
%REFERENCES:
% [1] "Ionospheric Correction Algorithm for Galileo Single-Frequency Users"
%      - European GNSS (Galileo) Open Service
%
%AUTHOR(s):
%Luigi De Maria, Matteo D'Addazio, 2022
%

%% Main Code

%Constants
DR = pi/180;                %Conversion Factor: deg->rad
RD = 180/pi;                %Conversion Factor: rad->deg
RE = 6371.2;                %Earth Mean Radius [km]

%Point Altitude
hs = sqrt(s^2 + perigee.rp^2) - RE;

%Trigo Values
tandels = s/perigee.rp;
cosdels = 1 / sqrt(1 + tandels^2);
sindels = tandels * cosdels;
%Latitude [rad]
sinlats = sin(perigee.latp*DR)*cosdels +...
          cos(perigee.latp*DR)*sindels*perigee.cossigp;
coslats = sqrt(1 - sinlats^2);
lats = atan2(sinlats, coslats) * RD;
%Longitude [rad]
sindiff = sindels * perigee.sinsigp * cos(perigee.latp*DR);
cosdiff = cosdels - sin(perigee.latp*DR)*sinlats;
lons = (atan2(sindiff, cosdiff) + perigee.lonp*DR) * RD;

end